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Abstract 

The well-known Hill’s averaging theorems for stresses and strains as well as the so-called 
Hill-Mandel principle of macrohomogeneity are essential ingredients for the coupling and the 
consistency between the micro and macro scales in multiscale finite element procedures (FE^). 
We show in this paper that these averaging relations hold exactly under standard finite element 
discretizations, even if the stress field is discontinuous across elements and the standard proofs 
based on the divergence theorem are no longer suitable. The discrete averaging results are 
derived for the three classical types of boundary conditions (affine displacement, periodic and 
uniform traction boundary conditions) using the properties of the shape functions and the 
weak form of the microscopic equilibrium equations. The analytical proofs are further verified 
numerically through a simple finite element simulation of an irregular representative volume 
element undergoing large deformations. Furthermore, the proofs are extended to include the 
effects of body forces and inertia, and the results are consistent with those in the smooth 
continuum setting. This work provides a solid foundation to apply Hill’s averaging relations in 
multiscale finite element methods without introducing an additional error in the scale transition 
due to the discretization. 


1 Introduction 

The vast majority of materials in nature as well as in engineering applications have underlying 
microstructures, and often, the length scale of the heterogeneites is much smaller than that of 
the system to be analyzed. In such cases, direct numerical simulations are typically prohibitive, 
and coarse-graining procedures have been developed to characterize the effective behavior of the 
material. A popular computational homogenization approach, that takes advantage of the sepa¬ 
ration of length scales, is the so-called FE^ method. This numerical strategy considers a (coarse) 
finite element discretization for the macroscopic domain, and evaluates the effective behavior of 
the material at each quadrature point through a representative volume element (RVE), where the 
microstructure is resolved. This approach is capable of dealing with general geometries, materials, 
and loading conditions, and takes into consideration the evolution of the microstructure. It has 
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called Hill’s theorems and Hill-Mandel principle of macrohomogeneity establish, in their original 
form, that the average strain and stress over the RVE are appropriate macroscopic quantities with 
which to describe the homogenized constitutive behavior, both in the linear and hnite kinematic 
setting (Zohdi and Wriggers, 2008; Hori and Nemat-Nasser, 1999). These average quantities depend 
exclusively on the value of the corresponding microscopic object at the boundary of the RVE; and 
they thus enable the formulation of a boundary value problem (with either Neumann, Dirichlet or 
periodic boundary conditions) from which the effective behavior may be obtained. 

To the best of the author’s knowledge, the above micro-macro relations and variational reformu¬ 


lations (Miehe et ah, 2002) and extensions to account for surfaces of discontinuities (Nemat-Nasser 


and Hori 

(2013), Section 2.4), or body forces 

or 

et al. 

2009 

Reina 

20111 Pham et al. 

2013 

d 


Ricker 


de Souza Neto et ah, 2015), all rely on the strong 


form of the equilibrium equation and successive application of the divergence theorem. However, 
these averaging relations are commonly used in finite element schemes, where the stress field may 
fail to be continuous, as is the case, for instance, for piecewise linear shape functions, and the 
balance equations are only satisfied in a weak sense. It is therefore, a priori, unclear whether Hill’s 
relations hold exactly, or only approximately, in a discrete setting. This is an important issue for 
error estimations in multiscale finite element methods. 

In this paper, the three fundamental averaging statements are shown to hold exactly for stan¬ 
dard finite element discretizations: (i) the volume-averaged deformation gradient relation; (ii) the 
volume-averaged stress relation; and (iii) the energy average relation or so-called Hill-Mandel prin¬ 
ciple of macrohomogeneity. The proofs are conducted initially in the static setting with no body 
forces, and then extended to the more general case where body forces and/or inertial effects are 
present. In each case, the three classical types of boundary conditions are considered: affine dis¬ 
placement, periodic, and uniform traction boundary conditions. The statements are derived from 
the properties of the shape functions and the weak form of the momentum balance equations, and 
the use of the divergence theorem is limited to continuous fields such as the deformation mapping. 
In the interest of generality, the proofs are conducted in the finite kinematic setting. Finally, a 
finite element simulation over a highly irregular RVE with a coarse mesh is employed to numerically 
verify the discrete averaging relations. Simple extension and simple shear deformation modes are 
considered using displacement, tractions or periodic boundary conditions. 


2 Notation 

The analyses are based on the Lagrangian formulation of continuum mechanics. In this setting, the 
deformation of a material point X in the reference configuration Qq C at time t is characterized 
by the deformation mapping x = (^(X, t). This mapping satishes the momentum balance equations. 
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which under sufficient smoothness, read 


( 1 ) 

( 2 ) 

( 3 ) 


V • P + B = 0, in Oq) 

(p = ip, on 9^0,1, 

PN = T, on 5f^o,2. 

where P is the first Piola-Kirchhoff stress tensor, (f and T are the prescribed deformation mapping 
and traction respectively, and N is the outward unit normal to the domain in the undeformed 
configuration. For static loading, B = Bq represents the body forces, whereas B also includes the 
inertial forces in the dynamic setting, i.e. B = Bq — Furthermore, the deformation gradient 
will be denoted by F = V(^, where V represents the material gradient. As usual, it is required 
that Sffo = U 0 ^ 0,25 and 9^0,1 H = 0- In order to distinguish between the micro- and 

macro-variables, the superscript M will be employed for the latter, whereas no superscript will be 
used for the microscopic quantities. 

The macroscopic fields will often result as the average of the corresponding microscopic objects. 
This average operation will be written as 

where |ffo| denotes the volume associated to Oq and dV is the volume differential. Where needed, 
the surface differential will be denoted as dS. 

Some derivations in the following sections will make use of standard indicial notation and 
Einstein summation convention. Lower case indices will then be used to refer to the deformed 
configuration and upper case indices for the reference configuration. 


3 Problem setting 

The multiscale finite element method FE^ solves a boundary value problem for an RVE at each 
quadrature point of the macroscopic scale. Different types of boundary conditions can be imposed 
at the RVE in order to couple the micro- and macro-solution: affine displacement, periodic or 
traction boundary conditions. 

The linear displacement boundary conditions read 

(^(X) = (p^ + F^X, on dQo, (5) 


where the macroscopic displacement field is often obviated in classical static analyses with 
divergence-free stresses, as a rigid body translation leaves the results unaltered, cf. [Nemat-Nasser 


and Hori (2013) Chapter 1. However, in the presence of body forces or inertial, it is both, physically 


and mathematically meaningful to inform the RVE of the macroscopic translations, the rotations 


already being included in the deformation gradient (Ricker et al., 2009 Reina 2011; de Souza Neto 


et al., 2015) 


Another boundary condition which is frequently employed is uniform traction boundary condi¬ 
tion, 

on cAIq. (6) 


T = P’^N, 
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It is well known that (Suquet 1987 Peric et ah, 2011) affine displacement boundary conditions 
result in stiffer solutions, whereas uniform traction boundary conditions are the most compliant. 
An intermediate behavior can be achieved with periodic boundary conditions 


ifiX) = if^ + F^X + if, 


on dQo, 


(7) 


where stands for the fluctuation field, which is periodic along each pair of parallel sides (faces in 
three dimensions). Similar to the case with affine boundary conditions, the macroscopic translation 
is unnecessary for static problems with no body forces. 


3.1 Finite element discretization and equilibrium equations 

Once a coupling strategy between the micro and macro scale is chosen, the multiscale FE^ problem, 
obeying Eqs. (§-(§ at both scales, is resolved via a hnite element discretization, not necessarily of 
the same type for the micro and macro problem. Here we consider conventional finite element 
discretizations ([Hughes 2012) for the RVE of the form 




( 8 ) 


where {a} represents the set of nodes, with associated degrees of freedom and Na are the 
corresponding shape functions, smooth within each element. These shape functions are required to 
have local support (each Na vanishes over any element not containing the node a) and to have the 
Kronecker delta property, Na(X.b) = ^ab- Furthermore, they shall satisfy the properties of partition 
of unity and linear field reproduction, i.e., 

J]iV,(X) = l, and 

a 

^iV,(X)X„j = Xj, 

a 

which allow an exact representation of arbitrary rigid body motions and uniform deformations. 
These requirements for the shape functions will be essential for the derivations in Sections and 


(9) 

( 10 ) 



• {c} boundary nodes 
. {b} interior nodes 
{a}={b}+{c} 


Figure 1: Schematic representation of a finite element discretization, where the nodes are colored 
according to their location, in the interior or at the boundary of the domain. 
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The weak form of the equilibrium equations associated the finite element discretization can then 
be readily obtained from the strong form, cf. Eq. Q. Towards that goal, it will result convenient to 
separate the node set {a} into interior nodes (body nodes) {6} and boundary nodes (at the contour 
of the domain) {c}, cf. Fig. and make use of the fact that the shape functions Nt, have a zero 
value at the boundary of the domain (by the local support property of the shape functions). 
Then, the weak form of the balance equations reads 




(V • P + B) • = 0 


E 


5ipa 




{PijNa,J - BiNa, 


dV — 5(pc 


PijNjNcdS = 0 . 


IdUo 


( 11 ) 

Equation ( [IT| ) shall be satisfied for any admissible variation of the nodal positions and, in 
particular, for the variation of any specific interior node b, while setting a null value for the variation 
of all the other nodes. Equivalently, 


/ {PijNb^j — BiNb) dV = 0, for all the interior nodes b. (12) 

J 

The variations of the nodes {c} depend on the boundary conditions used. For affine displacement 
boundary conditions, cf. Eq. Q, 6(p^ =0, and therefore no additional equations follow. For periodic 
boundary conditions, cf. Eq. Q, and the equilibrium equations for the boundary nodes 

read 

5ifci j {PijNc,j — BiNc) dV = 0 for all periodic. (13) 

c •'^0 

In the absence of inertial and body forces, this equation is equivalent to anti-periodic boundary 
surface traction. Finally, for uniform traction boundary conditions, cf. Eq. Q, the weak form of 
the balance law for the boundary nodes is 


[ {PijNcJ - BiN,) dV- [ P^NjN, dS = 0. (14) 

J r2Q J 

These equations may be readily simplified for the case with no body forces or inertia (B = 0), 
which is studied first. 

4 Discrete averaging results in the static case with no body forces 

The averaging statements for a representative volume element were initially developed for systems 
with negligible inertial and body forces and consist on the following three statements 

• Averaging theorem for the deformation gradient: for any F compatible (i.e. F = V(^) and 
(p = F'^X on Sflo, with ip periodic or null, the macroscopic deformation gradient F^ is 
equal to the average of its microscopic counterpart: F^ = (F). 

• Averaging theorem for the first Piola-Kirchhoff stress tensor: for any P in equilibrium (i.e. V • 
P = 0) and T = P^N on the macroscopic stress tensor P^ is equal to the average of 
the microscopic analog: P'^ = (P). 
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• Hill-Mandel principle: for any F compatible and any P in equilibrium, not necessarily related 
to each other, and for any of the three standard types of boundary conditions (affine dis¬ 
placements, periodic or stress boundary conditions), the following equality holds: (P : F) = 
(P) : (F). In view of the previous two relations, this principle establishes the energy rate 
consistency between the micro and macro scale. 


In this section, we show that these averaging results hold exactly under a finite element dis¬ 
cretization. This implies that only the nodes at the boundary are required to satisfy the displace¬ 
ment boundary conditions, and the equilibrium equations, both, in the body and with the external 
tractions, are only satisfied weakly, cf. Eqs. (12)-(14). Furthermore, for the Hill-Mandel principle, 
the deformation gradient F and the stress tensor P will not be required to be related to each other, 
but they shall result from an identical finite element discretization (same set of nodes and associated 
shape functions). The proofs are shown below for the three types of boundary conditions. 


4.1 Linear displacement bonndary conditions 

The first case considered is that of affine displacement boundary conditions of an RVE, in accordance 
with the macroscopic deformation gradient F^. In that case the boundary nodes {c} are required 
to satisfy 

99 ,, = f!^X,q. (15) 


Averaging statement for the deformation gradient 


Equation (15) is sufficient to show that the macroscopic deformation gradient tensor is the volume 
average of its microscopic analogue over the RVE. Since the displacement held is continuous over 
the domain and smooth in each element of the hnite element discretization, the divergence theorem 
can be directly applied to over Hq. Then, 


[ dV= [ if^Nj dS = y] ifai [ NaNj dS = y^ip,i [ . 

Jno JdUo n "''9^0 r "'aOo 


NcNj dS, 


(16) 


where the sum has been simplihed to the boundary nodes, since the shape functions associated to 
the interior nodes, N^, have a zero value at dfio- Next, we make use of the displacement boundary 
conditions on the node set {c}, cf. Eq. (15), and, again, use the fact that N/, has a zero value at 

dQo 


[ ipldV = F^Y.^<^Q [ N,NjdS = F^Y.^-Q [ 

Juq „ Jano n JdUo 


NaNjdS. 


(17) 


Finally, by the linear representation property of shape functions, cf. Eq. (10), and the application 
of the divergence theorem, the sought-after result is obtained 


[ dV = F^ [ XqNj dS = F^[ Xq^j dV = (18) 

or equivalently, 

F^ = {F) = j^[ FdV. (19) 

l^^ol Jno 

In analogy with the continuum strain averaging result, the above proof is purely kinematical in 
nature and it is independent of the microscopic equilibrium equations. 
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Hill-Mandel principle 


Next, we proceed to prove the so-called Hill-Mandel principle for a compatible rate of deformation 
gradient F = V(^^, and a stress tensor P in equilibrium, i.e. satisfying Eq. (12) with Bi = 0. The 
average microscopic (virtual) power can be written as 

[ Pij^ljdV=l Pijy^ipaiNa,jdV = y^ipbi [ P^JNt,JdV+[ Pijy^ifciN^jdV, {20) 
Qn nn ^ L iln 


where the sum over all nodes {a} has been divided into the sum over the interior nodes {b} and the 
boundary nodes {c}. The sum over {6} vanishes by the weak form of the equilibrium equations, 
cf. Eq. (12) with Bi = 0; and the sum over {c} can be rewritten, applying the boundary conditions 
given by Eq. (15), as 


[ Pij^ljdV=f PijY,Ff^X,QN^^jdV = [ P^jN,^jdV 


rpM 

^iQ- 


( 21 ) 


The sum can then be extended to the set of all nodes by the equilibrium equations of the inner 


nodes, cf. Eq. (12) with Bi = 0, from which it follows that 

[ PijY^XaQNajdV 
JUn „ 


P^JFiJdV = 


'Oo 


rpM 

^iQ- 


( 22 ) 


Additionally, from the properties of the shape functions, cf. Eq. (10), it is readily obtained that 

Y,^aQNa,J = XQ,j = 6Qj, (23) 


and therefore 


PijFijdV = 


' f2o 


PrjdV 


.JUq 




(24) 


By making use of the previously derived relation for the deformation gradient, i.e. = (Fij), we 
have 


(P : F) = (P) : (F). (25) 

This completes the proof of Hill-Mandel principle under a finite element discretization with linear 
displacement boundary condition. 

4.2 Periodic boundary condition 

The Dirichlet boundary conditions considered in the previous section may be relaxed, i.e. softer 
response, by allowing a periodic fluctuation field ip. Under these considerations, the boundary 
nodes are required to satisfy the following relationship, 

Pci = FIqXcQ + Pci, with Pci periodic. (26) 
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Averaging statement for the deformation gradient 


The proof of the average theorem for the deformation gradient follows in a very similar manner to the 


case with Dirichlet boundary conditions. Indeed, by making use of Eq. (16) and the displacement 


constraints on boundary nodes, cf. Eq. (26), it follows that 


[ vljdV = y^^cil N,NjdS = F^y^X,Q [ N,NjdS + y^^ci ! N.NjdS. ( 27 ) 
ioo c ddo-Q ^ Jano c ddfio 

The last term on the right hand side of the equation vanishes as each pair of boundary edges have an 
identical displacement fluctuation and opposite outward normals. Then, by extending the sum over 
{c} to the full set of nodes (zero value of Nf, on SOq), making use of the linear field reproduction 


property of the shape functions, cf. Eq. (10), and the divergence theorem, it follows that 


[ dV = F^Y. f dS = F!^ [ XqNj dS = F^ [ Xq^j dV = \Qo\Ff^. 

J rio d 5^0 'd d^o ■/ Qq 

(28) 

Thus, the classical averaging statement of the macroscopic deformation tensor is recovered, = 
(F), and is again independent of microscopic equilibrium equations. 


Hill-Mandel principle 

Next, we proceed to prove Hill’s energy consistency relation for a deformation rate F compatible 
and deriving from a periodic mapping and for a stress tensor P satisfying the weak equilibrium 


equations given by Eqs. (12) and (13) with Bi = 0. Similarly to Eq. (20), the average of the 


microscopic energy rate obtained from the finite element solution is given by 




PupljdV = 


b 


vfio 


PijNb,jdV 


^ib + 


/ PiJ VciNc, 


jdV, 


(29) 


where the set of nodes {a} has been divided into interior nodes {6} and nodes at the contour {c}. 
The term associated to the inner nodes vanished by the equilibrium equations. Then, by making use 


of the periodic boundary conditions, cf. Eq. (26), and the anti-periodicity of PijNc^j, cf. Eq. (13), 
one obtains 

[ PupIj dV= [ Pu V dV= [ Pu V (f^X,q + dV 

Juo Jno c dfio c ^ 


[ PuTx,QNc,jdv]F!^ + T^,, 

IJno r. J 


PijN.jdV 


.Jflo 


[ Puy^X,QN,^jdV 

.JUo r 


^iQ- 


(30) 


The proof then follows in an identical manner to the case with affine boundary conditions, see 


derivations after Eq. (21), leading as well to (P : F) = (P) : (F). This completes the proof of 
Hill-Mandel principle in discrete setting with periodic boundary conditions. 



















4.3 Uniform traction boundary conditions 

The third and last case considered is that of uniform traction boundary conditions, given by the 
macroscopic stress tensor 

T = P’^N. (31) 


Averaging statement for the first Piola-Kirchhoff stress tensor 

We proceed to show the volume-averaged relation of the stress tensor, under the boundary con¬ 
ditions given by Eq. (31). The volume integral of the microscopic stress over the RVE domain 
is 


[ Pij dV= [ P^qSqj dV= [ PiQ V XajNa,Q dV = y] X,j [ PiQN,,Q dV, 
J Qq ^0 ^0 n C ^0 


(32) 


where we have used the property of the shape functions given by Eq. (23) and the equilibrium 
equations for interior nodes {b}, cf. Eq. (12) with Bi = 0. Then, by substituting the weak form of 


NqN.dS. 


(33) 


the traction boundary condition, cf. Eq. (14) with Bi = 0, one obtains 

[ PudV = Y,Xcj [ P!^NQN,dS = I 

J rjo c 5r2o Q ^ ^ 

We may then extend the sum to all nodes {a}, as the shape functions associated to nodes {b} are 
zero at dflo- This results in 

[ Pu dV = P!^Yl f dS = P!^ [ NqXj = p!^ [ SjQ dV = P!j\no\ (34) 

Jflo a ddflo JdQo jflo 

where we have further used the exact linear representation property of the shape functions, cf. Eq. (10), 
and applied the divergence theorem. The well-known averaging theorem is then recovered exactly 


= (P) = 


1 


l^^ol 




PdV, 


(35) 


and follows from the equilibrium equations, in this case, in weak form. 


Hill-Mandel principle 

Einally, we derive the Hill-Mandel principle for the finite element problem with stress boundary 
conditions. Similar to the previous section, we make use of the governing equations corresponding 
to the node sets {6} and {c}, cf. Eqs. (12) and (0 with Bi = 0, respectively, and obtain 


(36) 


f Pij^lj dV= f Pij Y, ^aiNa,J dV = Y Pbi [ PijNb,j dV + Y 
4 f2o J Qo a ^ 4 Oo J flo 

= V / P!jNjN, dS = P!jY Vci [ NjN, dS. 

The shape functions associated to nodes {6} are zero at d^o and the last sum of the above equation 
may thus be extended to the full set of nodes {a}. It then follows that 

[ Pu^Ij dV = Pf^Y f dS = P!jY f dV = P!j [ dV, (37) 
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where the divergence theorem has been applied. By substituting the previously derived relation 
Pfj = {Pij), we obtain the sought-after result 

(P : F) = (P) : (F). (38) 


5 Discrete averaging results in the presence of body forces or in¬ 
ertia 

In this section, we extend the previous results to account for the presence of body forces or inertial 
terms, which are both considered in the term B. The boundary value problem for the RVE with 
periodic displacement boundary conditions and traction boundary conditions are considered sepa¬ 
rately, as well. The case with affine Dirichlet boundary conditions is omitted here, as its derivation 
is analogous to that of the periodic scenario. 

5.1 Periodic boundary condition 

In order to consider the work performed by the inertial and body forces, the boundary conditions 
on the microscopic domain read, for the parodic case, as 

Vci = + F^Xcq + ifci ■ (39) 

where informs the RVE of the macroscopic translation. 


Averaging statement for the deformation gradient 


As mentioned through the narrative, the averaging results for the deformation gradient are purely 
kinematical in nature and independent of the equilibrium equations. They thus hold as well in the 
presence of body forces or inertial terms. Note that the term has an identical treatment to ip 


in the derivations and therefore, the proof remains the same as that of Section 4.2 


Hill-Mandel principle 

A micro-macro relation analogous to the classical Hill-Mandel principle has been obtained for non- 


zero value of the body forces or inertia ( 

Molinari and Mercier 

2011 

Pham et al. 

2013 

de Souza Neto et al. 

2015 

). It reads 


2001 Ricker et al., 2009 Reina 


(P : F - B • (^) = (P - B ® X) : (F) - (B) • 


M 


(40) 


and of course, it reduces to the standard relation, (P : F) = (P) : (F), when the body forces and 
inertia terms vanish. 

In the following, we proceed to prove that this relation holds exactly under a finite element 
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discretization. The left hand side of the equation reads 

f (PiJ^ij - Bi^^) dV= [ Pu V ^aiNa,J dV- [ B,y] ipaiNa dV 
Jqq J Qq ^ J Uq ^ 

= J2^bi [ {PijNk^j - BiNk) dV +Y,^ci [ {PijN.^j - BiN,) dV 
iJ Uo J Lifio 

= V [ {PijNcJ - B,N,) dV 

^ L-'rin 


(41) 

where we have divided the full set of nodes into boundary nodes and interior nodes and have applied 
the equilibrium equation for the interior nodes, cf. Eq. (12). Next, we make use of the boundary 
conditions given by Eq. (39), and obtain 


[ [PijiPlj dV - B.if':) dV = y^ [ {P^JN^,J - BiN,) dV 

Jno ^ c 

+ \ [ (PaiVc, J - B,N,) dV 

+ Pci 


.J Uo 


.Jflo 

{PijNcJ - BiNc) dV 


(42) 


where the last term on the right hand side vanishes due to the periodicity of pci and anti-periodicity 
of the term in brackets, cf. Eq.(13). Then, by using Eq. (12) and extending the remaining sums to 
the full set of nodes {a}, one obtains 

[ (PuptjdV-B,p>l) dV 

Jno ^ J 


[ {PijNa,J - BiNa) dV 


[ {P^JNa,J - BiNa) dV 

_J np 

a 

.J rio 


= pr 


. 4 fig 




+ Fi 


. Jno 


PiJ Y. ^a,jXQa -B^Y NaXQa) dV 


Einally, the use of properties of the shape functions, cf. Eqs. and ([lol), delivers 


' Og 


(Pijp'lj - Bip'i) dV = [ {Pu-BiXj)dV 

^ J Uno 


pM 

d^iJ ~ 


BidV 


-Jno 


Pf 


(43) 


(44) 


where we have further used the fact that Eq. (§ implies X)a = 0- This above result, together 
with the already derived expression = (F), completes the desired proof. 


5.2 Uniform traction boundary condition 


When inertial or body forces are considered, a modified averaging statement for the the first Piola- 


Kirchhoff stress tensor has been obtained in the smooth case ( 

Molinari and Mercier 

2001 

Ricker 

et al. 

2009 

Reina 

2011 

Pham et al. 

2013 

de Souza Neto et ah, 2015). It reads 


pM ^ ^p_B0X). (45) 
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We show below that this equivalence holds exactly under a finite element discretization. 


Averaging statement for the first Piola-Kirchhoff stress tensor 


Using the properties of the shape functions, cf. Eqs. (10) and (23), and the eqnilibrium eqnations 
for interior nodes {6}, cf. Eq. (12), the right hand side of Eq. (45) reads 


' Qo 


{Pu-B,Xj)dV = / {P,q5qj - BiXj)dV = 


iQo 


[ f V XajNa,Q - BiV] NaXaj] dV 

\ a a / 


XaJ {PiQNa,Q - B^Na) dV = ^cj {PiQNc,Q - B,N,) dV. 

(46) 


'Oo 


' Qo 


Then, by substitnting the governing eqnations of boundary nodes {c}, cf. Eq. (14), we have 

^ {Pu - B,Xj) dV = P^Y. f 


' Oo 


'aoo 


(47) 


The shape functions associated to nodes {b} are zero at c?Uo and therefore the sum in the last term 
of the above eqnation may be extended over all nodes. Then, by the linear representation of the 
shape functions, cf. Eq. (10), and application of the divergence theorem, the sought-after result is 
obtained 


[ {Pu - B,Xj) dV = P!^y] XaJ [ NgNa dS = Pf^ [ NqXj dS = 
J n ^^0 


(48) 


Hill-Mandel principle 


In this section, we proceed to derive Eq. (40) for the finite element solution of an RVE with uniform 
traction bonndary condition. In this case, the finite element solution may be non-unique due to 
the undetermined value of the translations for the static problem, among other possible sources 
of non-uniqueness. Yet, these translations are important in the presence of body forces, as they 
contribute to the work these forces perform. It then resnlts convenient to separate the nodal 
deformation mapping Lpa f^e macroscopic translation, and the fluctuation around that valne. 


Eqnivalently, 


= 




which satisfies = X(p . Then, using the weak form of the 


equilibrium equations associated to nodes {6}, cf. Eq. (12), it follows that 



' Qo 


P^J Y, ^aiNa,J dV- ^aiNa + ) dV 


Y dV - B,Na) dV- [ Bapf dV 

a J Qo J Qo 


flo 

■M 


Y^^ci iPijNa,jdV-B,Na)dV-ipf B,dV. 




Oo 


(49) 
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Then, by substituting the governing equations of boundary nodes {c}, cf. Eq. (14), we obtain 




dV = Y^ [ pIqNqNJS - iff [ B, dV 
= NQN,dS-^f^[ B,dV. 


The shape functions associated to nodes {b} are zero at and, therefore 


'Oo 


(Pupt - BiPi) dV = P^Yl f - P 

^ 'a JdQo 


Bi dV. 


IQq 


Finally, application of divergence theorem gives, 

(Pupt - Bip’l) dV = P!j [ ^'iNj dS-^f [ Bi dV 

^ ' Idflo 


' Oo 


= P!jl ^IjdV-^, 

J r2o 


' Oo 
Bi dV. 


'Qo 


( 50 ) 


(51) 


(52) 


Substituting the previously obtained relation = (P — B 0 X), we obtain the desired relation 

(P : F - B • V?) = (P - B ® X) : (F) - (B) • ip^. (53) 


6 Numerical example 

In this section, the volume-averaged relations are verified numerically over an RVE for the three 
types of boundary conditions considered (affine, periodic, and traction boundary conditions), and 
two modes of deformation (simple extension and simple shear). Figure j^a) shows the domain of 
the RVE, which is a unit square composed of three materials: Al, Cu and Ni. These materials are 
arranged with no geometric symmetry, providing a general and irregular form of RVE. Furthermore, 
a coarse mesh is intentionally employed to highlight the validity of the results, even when far away 
from the smooth exact solution. In particular, the microscopic domain is divided into 8 identical 
linear triangular elements, as shown in Figurej^b), and linear shape functions are used for the finite 
element discretization. Simulations are performed using COMSOL multiphysics software, where 
hyperelastic constitutive relations (neo-Hookean) are used for the three materials, whose properties 
are listed in Table The details of the simulations are provided below along the text. 



(a) RVE model 


(b) Mesh distribution 


Figure 2: Geometry of RVE and undeformed mesh distribution. 
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Table 1: Materials properties for Al, Cu, and Ni from COMSOL material library. 



Al 

Cu 

Ni 

Density p (kg/m®) 

2700 

8960 

8900 

Lame constant A (GPa) 

60.49 

95.15 

136.4 

Lame constant p (GPa) 

25.93 

44.78 

80.15 


6.1 AfRne and periodic displacement boundary condition 

We begin by considering the plane-strain extension of the the RVE previously described according 
to the macroscopic deformation gradient F|f, F^] = [1.2, 0; 0,1]. Computationally, 

this macroscopic deformation is applied in 20 uniform steps, and at each of these steps, the stress 
and deformation gradient of each element is computed. Tables and show the values of each 
component of F and P for each element in the final configuration, represented in Figures j^a) and 
(b), for the affine and the periodic boundary conditions, respectively. 



(a) AfRne displacement b.c. 



(b) Periodic b.c. 


Figure 3: Deformed meshes for simple extension mode under affine and periodic boundary condi¬ 
tions (b.c.). 


As the results of Tables [^b) and ib) indicate, the average relation F^ = (F) is satisfied 
exactly, up to machine precision, for both affine and periodic boundary conditions. That is. 


1 

W\ 


f 

’|ELl(^ll)e lEl=liFl2)e 


1.2 0 

/ FdV = 


= 


J Uq 

iELl(^2l)e |ELl(^22)e_ 


0 1 


(54) 


where the sum is performed over the elements of the triangular mesh. A similar calculation can 
be carried out at each time step for the deformation gradient and the stress tensor, allowing us to 
verify the Hill-Mandel principle. For the linear displacement boundary conditions, the volumetric 
average of the strain energy stored during the deformation process is 


(P : F)dt 


'to 



( / P : Fdv) dt = 3.758512926F9 Pa, 

Vwol JOo / 


(55) 
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while the energy evaluated from the macroscopic (averaged) quantities leads to 

[ {P):{F)dt= [ {P):F^dt= [ ([ P : dt = 3.758512926E9 Pa, (56) 

Jto Jto Jto VI“o| Jno J 

which has, as expected, an identical value. Similarly, for the periodic boundary condition, the 
Hill-Mandel principle is verified exactly, with (P : F) = (P) : (F) = 3.638415107E9 Pa. 

A similar verification analysis is performed for a simple shear mode according to the macroscopic 
deformation gradient F'^ = [F^, Ff^', Eji > ; 0,1]. The corresponding deformed meshes 

of RVE are shown in Figure @ for the affine and periodic boundary conditions, respectively, 
and the components of the stress and deformation gradient tensor for the final configuration are 
recorded in Tables and for each element of the mesh. As indicated in the tables, the average 
deformation gradient coincides for both types of boundary conditions with the imposed macroscopic 
deformation, verifying the averaging relation F^ = (F). 




(a) Affine displacement b.c. 


(b) Periodic b.c. 


Figure 4: 
(b.c.). 


Deformed meshes for simple shear mode under affine and periodic boundary conditions 


Similarly, Hill-Mandel principle holds exactly under the finite element discretization. For linear 
displacement boundary condition, //^(P '■ F) dt = //^(P) : (F) dt = 9.569752719E8 Pa, and for 
periodic boundary conditions, //^(P :F) dt = //^(P) : (F) dt = 8.395522388E8 Pa. 

6.2 Uniform traction boundary conditions 

Next, we examine the averaging relations when the RVE is subjected to uniform traction boundary 
conditions. In particular, the following examples consider an applied macroscopic traction of P'^ = 
P^, P^] = [l, 0; 0, 0] GPa for a simple extension mode and P^=[P^Y) -^ 12 ! = 

1; 1, 0] GPa for a simple shear mode. The loading step size is set to O.lGPa when the non-zero 
stress component is smaller than 0.4 GPa, and to O.OlGPa in the following range. Figure [^a) and 
(b) shows the deformed meshes of RVE for both modes of deformation. 
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(a) Simple extension mode 


(b) Simple shear mode 


Figure 5: Deformed meshes for simple extension mode and simple shear mode under uniform 
traction boundary condition. 


The components of the stress tensor and the deformation tensor at the last step of deformation 
for all the elements have been listed in Tables and for simple extension mode and simple shear 
mode, respectively. The results indicate that the volume average of the stress tensor is exactly equal 
to P^, up to machine precision, verifying the discrete stress averaging result derived in Section 
Similarly, the Hill-Mandel principle is verihed numerically. For simple extension mode, the micro- 
and macro-work coincide and are equal to //^(P '■ F) dt = //^(P) : (F) dt = 3.726695734E6 Pa, 
whereas for the simple shear mode, //^(P :F) dt = f^^{F} ■ (F) dt = 1.192942091E7Pa. 

6.3 Discussion 

The numerical examples above confirm the validity of the stress and strain averaging relations as 
well as Hill-Mandel principle for a finite element discretization regardless of the size of the RVE 
or the mesh size. The energetic consistency relation, i.e. (P : F) = (P) : (F), is demonstrated for 
stresses and strains that belong to the same boundary value problem, as this is often its primary use. 
Yet, this relation only requires a compatible deformation gradient and a stress tensor in equilibrium 
that share the same finite element discretization. However, they are not required to be related to 
each other, as shown in the analytical derivations. 

7 Conclusion 

In summary, the averaging relations needed to establish the consistency between the two levels of 
description in computational homogenization methods (FE^) are shown to be exact under standard 
finite element discretizations. These are the stress and strain averaging relations as well as Hill- 
Mandel principle of energetic consistency. The proofs are performed in the finite kinematic setting, 
for generality, and are shown to hold for the three standard type of boundary conditions (affine 
displacement, uniform stress and periodic boundary conditions) for static or dynamic loading. The 
analytical proofs are further supported by numerical examples of a highly asymmetric RVE with 
a coarse finite element discretization. This work provides a solid foundation to the applications 
of Hill’s averaging theorems in RVE-based multiscale finite element methods with conventional 
Dirichlet, Neumann or periodic boundary conditions. 


4.3 
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Appendix A 


Table 2a: Components of stress tensor for simple extension mode under linear displacement bound¬ 
ary condition. 


Number 

Pii(Pa) 

Pi 2 (Pa) 

P 2 i(Pa) 

P22(Pa) 

1 

5.540445662129311E10 

-3.697775921110617E9 

-3.697775921110617E9 

1.872085738015028E10 

2 

4.564975169177103E10 

-2.311265268571417E9 

-2.311265268571417E9 

2.098423482999226E10 

3 

3.116295439474504E10 

7.745539252061372E8 

7.745539252061372E8 

1.996623282454363E10 

4 

3.694029850600000E10 

0 

0 

1.902985074600000E10 

5 

3.694029850600000E10 

0 

0 

1.902985074600000E10 

6 

4.271764261725496E10 

-7.745539252061372E8 

-7.745539252061372E8 

1.809346866745638E10 

7 

2.765210896757578E10 

7.475944449735475E8 

7.475944449735475E8 

1.488959713525095E10 

8 

2.421352284198325E10 

1.196070730162412E9 

1.196070730162412E9 

1.533834135993035E10 


Table 2b; Components of deformation gradient for simple extension mode under linear displacement 
boundary condition. 


Number 

Ell 

El 2 

E 21 

E 22 

1 

1.2 

-0.04613415672758626 

0 

0.9711642142698373 

2 

1.153865843272414 

0 

-0.028835785730162634 

1 

3 

1.153865843272414 

0.04613415672758626 

-0.028835785730162634 

1.028835785730163 

4 

1.2 

0 

0 

1 

5 

1.2 

0 

0 

1 

6 

1.246134156727586 

-0.04613415672758626 

0.028835785730162634 

0.9711642142698373 

7 

1.246134156727586 

0 

0.028835785730162634 

1 

8 

1.2 

0.04613415672758626 

0 

1.028835785730163 

Average 

1.200000000000000 

0 

0 

1 
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Table 3a: Components of stress tensor for simple extension mode under periodic boundary condi¬ 
tion. 


Number Pii(Pa) Pi 2 (Pa) P 2 i(Pa) P 22 (Pa) 


1 

2 

3 

4 

5 

6 

7 

8 


4.272877655616882E10 

4.272877655616881E10 

3.003952559993056E10 

3.003952559993057E10 

4.384107141206944E10 

4.384107141206944E10 

2.892723074402995E10 

2.892723074402994E10 


6.532528423931364E-7 

2.780854427286705E-7 

1.398135338753770E-6 

1.297020385082037E-6 

-5.423341730090870E-8 

-1.242786967781129E-6 

-8.994862469441367E-7 

-3.140181817261574E-8 


6.532528423931364E-7 

2.780854427286705E-7 

1.398135338753770E-6 

1.297020385082037E-6 

-5.423341730090870E-8 

-1.242786967781129E-6 

-8.994862469441367E-7 

-3.140181817261574E-8 


1.517259313986339E10 

1.517259313986339E10 

1.806628499750840E10 

1.806628499750839E10 

1.999341649449161E10 

1.999341649449161E10 

1.709972463684661E10 

1.709972463684661E10 


Table 3b: Components of deformation gradient for simple extension mode under periodic boundary 
condition. 


Number 

Pii 

Pi 2 

P 21 

P 22 

1 

1.152799641638063 

0 

0 

0.9809015450978975 

2 

1.152799641638063 

0 

0 

0.9809015450978975 

3 

1.152799641638063 

0 

0 

1.019098454902103 

4 

1.152799641638063 

0 

0 

1.019098454902103 

5 

1.247200358361937 

0 

0 

0.9809015450978975 

6 

1.247200358361937 

0 

0 

0.9809015450978975 

7 

1.247200358361937 

0 

0 

1.019098454902103 

8 

1.247200358361937 

0 

0 

1.019098454902103 

Average 

1.200000000000000 

0 

0 

1 


Table 4a: Components of stress tensor for simple shear mode under linear displacement boundary 
condition. 


Number 

Pii(Pa) 

Pi 2 (Pa) 

P 2 i(Pa) 

P 22 (Pa) 

1 

-2.130595379666259E9 

1.477835058908113E10 

1.477835058908113E10 

-4.634962905504009E9 

2 

-4.634962905504009E9 

1.477835058908113E10 

1.477835058908113E10 

-2.130595379666260E9 

3 

-1.399028340492011E9 

8.955223880000000E9 

8.955223880000000E9 

1.399028340492011E9 

4 

0 

8.955223880000000E9 

8.955223880000000E9 

0 

5 

0 

8.955223880000000E9 

8.955223880000000E9 

0 

6 

1.399028340492011E9 

8.955223880000000E9 

8.955223880000000E9 

-1.399028340492011E9 

7 

1.755118475952594E9 

5.590212526640496E9 

5.590212526640496E9 

9.450637946716006E8 

8 

9.450637946716006E8 

5.590212526640497E9 

5.590212526640497E9 

1.755118475952594E9 
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Table 4b; Components of deformation gradient for simple shear mode under linear displacement 
boundary condition. 


Number 

Pir 

Pi 2 

P 21 

P 22 

1 

1 

0.1843775168634644 

0 

0.9843775168634644 

2 

0.9843775168634644 

0.2 

-0.01562248313653562 

1 

3 

0.9843775168634644 

0.2156224831365357 

-0.01562248313653562 

1.015622483136536 

4 

1 

0.2 

0 

1 

5 

1 

0.2 

0 

1 

6 

1.015622483136536 

0.1843775168634644 

0.01562248313653562 

0.9843775168634644 

7 

1.015622483136536 

0.2 

0.01562248313653562 

1 

8 

1 

0.2156224831365357 

0 

1.015622483136536 

Average 

1 

0.2000000000000000 

0 

1 


Table 5a: Components of stress tensor for simple shear mode under periodic boundary condition. 


Number Pii(Pa) Pi 2 (Pa) P 2 i(Pa) P 22 (Pa) 


1 

2 

3 

4 

5 

6 

7 

8 


1.892653704593794E-6 

7.064054015297273E-7 

-6.602305766684191E-7 

4.633653368252543E-6 

4.928428739833904E-7 

3.883709277785474E-8 

-1.092893533105891E-6 

-3.537979237633326E-6 


7.835820896554924E9 

7.835820896554925E9 

8.955223880000004E9 

8.955223880000004E9 

8.955223880000000E9 

8.955223880000000E9 

7.835820896554928E9 

7.835820896554928E9 


7.835820896554924E9 

7.835820896554925E9 

8.955223880000004E9 

8.955223880000004E9 

8.955223880000000E9 

8.955223880000000E9 

7.835820896554928E9 

7.835820896554928E9 


4.117337246423158E-6 

1.536736099019807E-6 

-1.281624060558984E-6 

1.684227316024411E-6 

9.566949906492653E-7 

1.903017544449548E-6 

-1.001675136507230E-6 

-2.458595741300362E-6 


Table 5b: Components of deformation gradient for simple shear mode under periodic boundary 
condition. 


Number 

Pii 

Pi 2 

P 21 

P 22 

1 

1 

0.1488805970187594 

-0.05111940298124053 

1 

2 

1 

0.1488805970187594 

-0.05111940298124052 

1 

3 

1 

0.2511194029812406 

-0.05111940298124052 

1 

4 

1 

0.2511194029812406 

-0.05111940298124052 

1 

5 

1 

0.1488805970187594 

0.05111940298124053 

1 

6 

1 

0.1488805970187594 

0.05111940298124052 

1 

7 

1 

0.2511194029812406 

0.05111940298124052 

1 

8 

1 

0.2511194029812406 

0.05111940298124052 

1 

Average 

1 

0.2000000000000000 

0 

1 
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Table 6a: Components of stress tensor for simple extension mode under uniform traction boundary 
condition. 


Number Pii(GPa) Pi 2 (GPa) P 2 i(GPa) P 22 (GPa) 

1 1.039825026165646 0.04820528954458368 0.04866128855931545 -0.2416297163702212 

2 1.320203739079774 -0.03982502616563410 -0.03874395172276655 -0.04866128855930084 

3 0.6802710990595731 0.02932516104756559 0.02946016137955194 -0.1912627316990225 

4 0.9597001356950392 -0.04029986430498470 -0.03937749821610080 -0.03937749821606967 

5 1.034848686668402 0.03484868666840837 0.03433176606052552 0.03433176606052036 

6 1.300707598404531 -0.04322895004735811 -0.04357035273468100 0.2559592388690298 

7 0.7026611918618899 0.03821747693482153 0.03716704160995366 0.02792845493580790 

8 0.9617825230651917 -0.02724277367740260 -0.02792845493579822 0.2027117749793306 

Average 1.000000000000006 -0.00000000000000004 0.000000000000000001 0.0000000000000093 


Table 6b: Components of deformation gradient for simple extension mode under uniform traction 
boundary condition. 


Number Pii P 12 P 21 P 22 

1 1.004937304707261 -2.689342182920912E-4 8.747038856957929E-4 0.9969235353724314 

2 1.005765870068560 -1.097499579591548E-3 6.054301474451633E-4 0.9971928091106821 

3 1.005765870068560 5.280152905005938E-5 6.054301474451633E-4 0.9960071216113445 

4 1.007269717597589 -1.451045999979014E-3 5.567276539679379E-4 0.9960558241048217 

5 1.007549975212814 8.729795948439010E-4 -9.583423724183159E-5 0.9963105708813713 

6 1.008691889025950 -2.689342182920912E-4 -7.087987283018901E-4 0.9969235353724314 

7 1.008691889025950 2.172935797205065E-3 -7.087987283018901E-4 0.9955870601989383 

8 1.010812023294105 5.280152905005938E-5 -1.128860140708145E-3 0.9960071216113445 


Table 7a: Components of stress tensor for simple shear mode under uniform traction boundary 
condition. 


Number Pii(GPa) Pi 2 (GPa) P 2 i(GPa) P 22 (GPa) 

1 0.04560639290814648 0.9543878256031890 0.9543936070918476 -0.015047074289714631 

2 -0.01504707428965104 0.9543936070918378 0.9543878256031800 0.045606392908200080 

3 -0.06074680125388290 1.060613153272475 1.061397992038798 -0.068434256720068380 

4 0.03018748263543562 1.030187482635433 1.029820575266174 0.029820575266171635 

5 0.02982057526613236 1.029820575266177 1.030187482635437 0.030187482635397278 

6 -0.06843425672007669 1.061397992038797 1.060613153272471 -0.060746801253881974 

7 -6.77692342451637E-3 0.9546093951216100 0.9545899689704889 0.045390604878376570 

8 0.04539060487842260 0.9545899689704804 0.9546093951216026 -6.7769234244772734E-3 

Average 0.0000000000000013 1.000000000000000 1.000000000000000 0.0000000000000004 


20 



Table 7b: Components of deformation gradient for simple shear mode under uniform traction 
boundary condition. 


Number 

Ell 

El 2 

E 21 

E 22 

1 

1.000242207581661 

5.764948592700328E-3 

6.143346763082536E-3 

0.9998638094113029 

2 

0.9998638094113029 

6.143346763057936E-3 

5.764948592724816E-3 

1.000242207581661 

3 

0.9998638094113029 

0.01791376094498578 

5.764948592724816E-3 

0.9997780280404228 

4 

1.000156756804026 

0.01762081355226269 

5.390318334383442E-3 

1.000152658298764 

5 

1.000152658298764 

5.390318334358900E-3 

0.01762081355228730 

1.000156756804026 

6 

0.9997780280404227 

5.764948592700328E-3 

0.01791376094501020 

0.9998638094113029 

7 

0.9997780280404227 

0.01892022363786956 

0.01791376094501020 

1.000784490733306 

8 

1.000784490733307 

0.01791376094498578 

0.01892022363789369 

0.9997780280404228 
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